// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 - 2020 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.

#pragma once
#ifndef KERNEL_CUBATURE_SHUNN_HAM_DRIVER_HPP
#define KERNEL_CUBATURE_SHUNN_HAM_DRIVER_HPP 1

// includes, FEAT
#include <kernel/cubature/symmetric_simplex_driver.hpp>
#include <kernel/util/meta_math.hpp>

namespace FEAT
{
  namespace Cubature
  {
    /**
     * \brief Shunn-Ham driver class template
     *
     * This driver implements the Shunn-Ham rules for tetrahedra.
     * \see Lee Shunn, Frank Ham: Symmetric quadrature rules for tetrahedra based on cubic
     *      close-packed lattice arrangement; Journal of Compututational and Applied Mathematics,
     *      Volume 236 (2012), pp. 4348 -- 4364
     *
     * \tparam Shape_
     * The shape type of the element.
     *
     * \author Peter Zajac
     */
    template<typename Shape_>
    class ShunnHamDriver DOXY({});

    /// actual implementation
    template<>
    class ShunnHamDriver<Shape::Simplex<3>> :
      public DriverBase<Shape::Simplex<3>>
    {
    public:
      /// this rule is variadic
      static constexpr bool variadic = true;
      static constexpr int min_points = 2;
      static constexpr int max_points = 6;

      /// Returns the name of the cubature rule.
      static String name()
      {
        return "shunn-ham";
      }

      static int count(int param)
      {
        switch(param)
        {
        case 2:
          return 4;
        case 3:
          return 10;
        case 4:
          return 20;
        case 5:
          return 35;
        case 6:
          return 56;
        default:
          return 0;
        }
      }

      // note: the coords and weights are only given in double precision,
      // so we use 'double' for the parameters here...
      template<
        typename Weight_,
        typename Coord_,
        typename Point_>
      static void set(Rule<Shape::Simplex<3>, Weight_, Coord_, Point_>& rule,
        int k, double w, double x, double y, double z)
      {
        rule.get_weight(k) = Weight_(w);
        rule.get_coord(k, 0) = Coord_(x);
        rule.get_coord(k, 1) = Coord_(y);
        rule.get_coord(k, 2) = Coord_(z);
      }

      /**
       * \brief Fills the cubature rule structure.
       *
       * \param[in,out] rule
       * The cubature rule to be filled.
       */
      template<
        typename Weight_,
        typename Coord_,
        typename Point_>
      static void fill(Rule<Shape::Simplex<3>, Weight_, Coord_, Point_>& rule, int num_points)
      {
        switch(num_points)
        {
        case 2:
          set(rule,  0, 0.0416666666666667, 0.1381966011250110, 0.1381966011250110, 0.1381966011250110);
          set(rule,  1, 0.0416666666666667, 0.5854101966249680, 0.1381966011250110, 0.1381966011250110);
          set(rule,  2, 0.0416666666666667, 0.1381966011250110, 0.5854101966249680, 0.1381966011250110);
          set(rule,  3, 0.0416666666666667, 0.1381966011250110, 0.1381966011250110, 0.5854101966249680);
          break;

        case 3:
          set(rule,  0, 0.0079388558072015, 0.0738349017262234, 0.0738349017262234, 0.0738349017262234);
          set(rule,  1, 0.0079388558072015, 0.7784952948213300, 0.0738349017262234, 0.0738349017262234);
          set(rule,  2, 0.0079388558072015, 0.0738349017262234, 0.7784952948213300, 0.0738349017262234);
          set(rule,  3, 0.0079388558072015, 0.0738349017262234, 0.0738349017262234, 0.7784952948213300);
          set(rule,  4, 0.0224852072396435, 0.4062443438840510, 0.0937556561159491, 0.0937556561159491);
          set(rule,  5, 0.0224852072396435, 0.0937556561159491, 0.4062443438840510, 0.0937556561159491);
          set(rule,  6, 0.0224852072396435, 0.0937556561159491, 0.0937556561159491, 0.4062443438840510);
          set(rule,  7, 0.0224852072396435, 0.4062443438840510, 0.4062443438840510, 0.0937556561159491);
          set(rule,  8, 0.0224852072396435, 0.4062443438840510, 0.0937556561159491, 0.4062443438840510);
          set(rule,  9, 0.0224852072396435, 0.0937556561159491, 0.4062443438840510, 0.4062443438840510);
          break;

        case 4:
          set(rule,  0, 0.0011778457990783, 0.0323525947272439, 0.0323525947272439, 0.0323525947272439);
          set(rule,  1, 0.0011778457990783, 0.9029422158182679, 0.0323525947272439, 0.0323525947272439);
          set(rule,  2, 0.0011778457990783, 0.0323525947272439, 0.9029422158182679, 0.0323525947272439);
          set(rule,  3, 0.0011778457990783, 0.0323525947272439, 0.0323525947272439, 0.9029422158182679);
          set(rule,  4, 0.0078331114953146, 0.6165965330619370, 0.0603604415251421, 0.0603604415251421);
          set(rule,  5, 0.0078331114953146, 0.2626825838877790, 0.0603604415251421, 0.0603604415251421);
          set(rule,  6, 0.0078331114953146, 0.0603604415251421, 0.6165965330619370, 0.0603604415251421);
          set(rule,  7, 0.0078331114953146, 0.0603604415251421, 0.2626825838877790, 0.0603604415251421);
          set(rule,  8, 0.0078331114953146, 0.0603604415251421, 0.0603604415251421, 0.6165965330619370);
          set(rule,  9, 0.0078331114953146, 0.0603604415251421, 0.0603604415251421, 0.2626825838877790);
          set(rule, 10, 0.0078331114953146, 0.2626825838877790, 0.6165965330619370, 0.0603604415251421);
          set(rule, 11, 0.0078331114953146, 0.6165965330619370, 0.2626825838877790, 0.0603604415251421);
          set(rule, 12, 0.0078331114953146, 0.2626825838877790, 0.0603604415251421, 0.6165965330619370);
          set(rule, 13, 0.0078331114953146, 0.6165965330619370, 0.0603604415251421, 0.2626825838877790);
          set(rule, 14, 0.0078331114953146, 0.0603604415251421, 0.2626825838877790, 0.6165965330619370);
          set(rule, 15, 0.0078331114953146, 0.0603604415251421, 0.6165965330619370, 0.2626825838877790);
          set(rule, 16, 0.0169894863816447, 0.3097693042728620, 0.3097693042728620, 0.0706920871814129);
          set(rule, 17, 0.0169894863816447, 0.3097693042728620, 0.0706920871814129, 0.3097693042728620);
          set(rule, 18, 0.0169894863816447, 0.0706920871814129, 0.3097693042728620, 0.3097693042728620);
          set(rule, 19, 0.0169894863816447, 0.3097693042728620, 0.3097693042728620, 0.3097693042728620);
          break;

        case 5:
          set(rule,  0, 0.0003650077327565, 0.0267367755543735, 0.0267367755543735, 0.0267367755543735);
          set(rule,  1, 0.0003650077327565, 0.9197896733368800, 0.0267367755543735, 0.0267367755543735);
          set(rule,  2, 0.0003650077327565, 0.0267367755543735, 0.9197896733368800, 0.0267367755543735);
          set(rule,  3, 0.0003650077327565, 0.0267367755543735, 0.0267367755543735, 0.9197896733368800);
          set(rule,  4, 0.0023899278362944, 0.7477598884818090, 0.0391022406356488, 0.0391022406356488);
          set(rule,  5, 0.0023899278362944, 0.1740356302468940, 0.0391022406356488, 0.0391022406356488);
          set(rule,  6, 0.0023899278362944, 0.0391022406356488, 0.7477598884818090, 0.0391022406356488);
          set(rule,  7, 0.0023899278362944, 0.0391022406356488, 0.1740356302468940, 0.0391022406356488);
          set(rule,  8, 0.0023899278362944, 0.0391022406356488, 0.0391022406356488, 0.7477598884818090);
          set(rule,  9, 0.0023899278362944, 0.0391022406356488, 0.0391022406356488, 0.1740356302468940);
          set(rule, 10, 0.0023899278362944, 0.1740356302468940, 0.7477598884818090, 0.0391022406356488);
          set(rule, 11, 0.0023899278362944, 0.7477598884818090, 0.1740356302468940, 0.0391022406356488);
          set(rule, 12, 0.0023899278362944, 0.1740356302468940, 0.0391022406356488, 0.7477598884818090);
          set(rule, 13, 0.0023899278362944, 0.7477598884818090, 0.0391022406356488, 0.1740356302468940);
          set(rule, 14, 0.0023899278362944, 0.0391022406356488, 0.1740356302468940, 0.7477598884818090);
          set(rule, 15, 0.0023899278362944, 0.0391022406356488, 0.7477598884818090, 0.1740356302468940);
          set(rule, 16, 0.0041717565947791, 0.4547545999844830, 0.0452454000155172, 0.0452454000155172);
          set(rule, 17, 0.0041717565947791, 0.0452454000155172, 0.4547545999844830, 0.0452454000155172);
          set(rule, 18, 0.0041717565947791, 0.0452454000155172, 0.0452454000155172, 0.4547545999844830);
          set(rule, 19, 0.0041717565947791, 0.4547545999844830, 0.4547545999844830, 0.0452454000155172);
          set(rule, 20, 0.0041717565947791, 0.4547545999844830, 0.0452454000155172, 0.4547545999844830);
          set(rule, 21, 0.0041717565947791, 0.0452454000155172, 0.4547545999844830, 0.4547545999844830);
          set(rule, 22, 0.0079973222176259, 0.2232010379623150, 0.2232010379623150, 0.0504792790607720);
          set(rule, 23, 0.0079973222176259, 0.5031186450145980, 0.2232010379623150, 0.0504792790607720);
          set(rule, 24, 0.0079973222176259, 0.2232010379623150, 0.5031186450145980, 0.0504792790607720);
          set(rule, 25, 0.0079973222176259, 0.2232010379623150, 0.0504792790607720, 0.2232010379623150);
          set(rule, 26, 0.0079973222176259, 0.5031186450145980, 0.0504792790607720, 0.2232010379623150);
          set(rule, 27, 0.0079973222176259, 0.2232010379623150, 0.0504792790607720, 0.5031186450145980);
          set(rule, 28, 0.0079973222176259, 0.0504792790607720, 0.2232010379623150, 0.2232010379623150);
          set(rule, 29, 0.0079973222176259, 0.0504792790607720, 0.5031186450145980, 0.2232010379623150);
          set(rule, 30, 0.0079973222176259, 0.0504792790607720, 0.2232010379623150, 0.5031186450145980);
          set(rule, 31, 0.0079973222176259, 0.5031186450145980, 0.2232010379623150, 0.2232010379623150);
          set(rule, 32, 0.0079973222176259, 0.2232010379623150, 0.5031186450145980, 0.2232010379623150);
          set(rule, 33, 0.0079973222176259, 0.2232010379623150, 0.2232010379623150, 0.5031186450145980);
          set(rule, 34, 0.0155290955199223, 0.2500000000000000, 0.2500000000000000, 0.2500000000000000);
          break;

        case 6:
          set(rule,  0, 0.0001728852056023, 0.0149520651530592, 0.0149520651530592, 0.0149520651530592);
          set(rule,  1, 0.0001728852056023, 0.9551438045408220, 0.0149520651530592, 0.0149520651530592);
          set(rule,  2, 0.0001728852056023, 0.0149520651530592, 0.9551438045408220, 0.0149520651530592);
          set(rule,  3, 0.0001728852056023, 0.0149520651530592, 0.0149520651530592, 0.9551438045408220);
          set(rule,  4, 0.0016002774233247, 0.1518319491659370, 0.0340960211962615, 0.0340960211962615);
          set(rule,  5, 0.0016002774233247, 0.7799760084415400, 0.0340960211962615, 0.0340960211962615);
          set(rule,  6, 0.0016002774233247, 0.0340960211962615, 0.1518319491659370, 0.0340960211962615);
          set(rule,  7, 0.0016002774233247, 0.0340960211962615, 0.7799760084415400, 0.0340960211962615);
          set(rule,  8, 0.0016002774233247, 0.0340960211962615, 0.0340960211962615, 0.1518319491659370);
          set(rule,  9, 0.0016002774233247, 0.0340960211962615, 0.0340960211962615, 0.7799760084415400);
          set(rule, 10, 0.0016002774233247, 0.7799760084415400, 0.1518319491659370, 0.0340960211962615);
          set(rule, 11, 0.0016002774233247, 0.1518319491659370, 0.7799760084415400, 0.0340960211962615);
          set(rule, 12, 0.0016002774233247, 0.7799760084415400, 0.0340960211962615, 0.1518319491659370);
          set(rule, 13, 0.0016002774233247, 0.1518319491659370, 0.0340960211962615, 0.7799760084415400);
          set(rule, 14, 0.0016002774233247, 0.0340960211962615, 0.7799760084415400, 0.1518319491659370);
          set(rule, 15, 0.0016002774233247, 0.0340960211962615, 0.1518319491659370, 0.7799760084415400);
          set(rule, 16, 0.0027415662799705, 0.5526556431060170, 0.0462051504150017, 0.0462051504150017);
          set(rule, 17, 0.0027415662799705, 0.3549340560639790, 0.0462051504150017, 0.0462051504150017);
          set(rule, 18, 0.0027415662799705, 0.0462051504150017, 0.5526556431060170, 0.0462051504150017);
          set(rule, 19, 0.0027415662799705, 0.0462051504150017, 0.3549340560639790, 0.0462051504150017);
          set(rule, 20, 0.0027415662799705, 0.0462051504150017, 0.0462051504150017, 0.5526556431060170);
          set(rule, 21, 0.0027415662799705, 0.0462051504150017, 0.0462051504150017, 0.3549340560639790);
          set(rule, 22, 0.0027415662799705, 0.3549340560639790, 0.5526556431060170, 0.0462051504150017);
          set(rule, 23, 0.0027415662799705, 0.5526556431060170, 0.3549340560639790, 0.0462051504150017);
          set(rule, 24, 0.0027415662799705, 0.3549340560639790, 0.0462051504150017, 0.5526556431060170);
          set(rule, 25, 0.0027415662799705, 0.5526556431060170, 0.0462051504150017, 0.3549340560639790);
          set(rule, 26, 0.0027415662799705, 0.0462051504150017, 0.3549340560639790, 0.5526556431060170);
          set(rule, 27, 0.0027415662799705, 0.0462051504150017, 0.5526556431060170, 0.3549340560639790);
          set(rule, 28, 0.0025624627752218, 0.2281904610687610, 0.2281904610687610, 0.0055147549744775);
          set(rule, 29, 0.0025624627752218, 0.5381043228880020, 0.2281904610687610, 0.0055147549744775);
          set(rule, 30, 0.0025624627752218, 0.2281904610687610, 0.5381043228880020, 0.0055147549744775);
          set(rule, 31, 0.0025624627752218, 0.2281904610687610, 0.0055147549744775, 0.2281904610687610);
          set(rule, 32, 0.0025624627752218, 0.5381043228880020, 0.0055147549744775, 0.2281904610687610);
          set(rule, 33, 0.0025624627752218, 0.2281904610687610, 0.0055147549744775, 0.5381043228880020);
          set(rule, 34, 0.0025624627752218, 0.0055147549744775, 0.2281904610687610, 0.2281904610687610);
          set(rule, 35, 0.0025624627752218, 0.0055147549744775, 0.5381043228880020, 0.2281904610687610);
          set(rule, 36, 0.0025624627752218, 0.0055147549744775, 0.2281904610687610, 0.5381043228880020);
          set(rule, 37, 0.0025624627752218, 0.5381043228880020, 0.2281904610687610, 0.2281904610687610);
          set(rule, 38, 0.0025624627752218, 0.2281904610687610, 0.5381043228880020, 0.2281904610687610);
          set(rule, 39, 0.0025624627752218, 0.2281904610687610, 0.2281904610687610, 0.5381043228880020);
          set(rule, 40, 0.0048920019729205, 0.3523052600879940, 0.3523052600879940, 0.0992057202494530);
          set(rule, 41, 0.0048920019729205, 0.1961837595745600, 0.3523052600879940, 0.0992057202494530);
          set(rule, 42, 0.0048920019729205, 0.3523052600879940, 0.1961837595745600, 0.0992057202494530);
          set(rule, 43, 0.0048920019729205, 0.3523052600879940, 0.0992057202494530, 0.3523052600879940);
          set(rule, 44, 0.0048920019729205, 0.1961837595745600, 0.0992057202494530, 0.3523052600879940);
          set(rule, 45, 0.0048920019729205, 0.3523052600879940, 0.0992057202494530, 0.1961837595745600);
          set(rule, 46, 0.0048920019729205, 0.0992057202494530, 0.3523052600879940, 0.3523052600879940);
          set(rule, 47, 0.0048920019729205, 0.0992057202494530, 0.1961837595745600, 0.3523052600879940);
          set(rule, 48, 0.0048920019729205, 0.0992057202494530, 0.3523052600879940, 0.1961837595745600);
          set(rule, 49, 0.0048920019729205, 0.1961837595745600, 0.3523052600879940, 0.3523052600879940);
          set(rule, 50, 0.0048920019729205, 0.3523052600879940, 0.1961837595745600, 0.3523052600879940);
          set(rule, 51, 0.0048920019729205, 0.3523052600879940, 0.3523052600879940, 0.1961837595745600);
          set(rule, 52, 0.0061048561067518, 0.1344783347929940, 0.1344783347929940, 0.1344783347929940);
          set(rule, 53, 0.0061048561067518, 0.5965649956210169, 0.1344783347929940, 0.1344783347929940);
          set(rule, 54, 0.0061048561067518, 0.1344783347929940, 0.5965649956210169, 0.1344783347929940);
          set(rule, 55, 0.0061048561067518, 0.1344783347929940, 0.1344783347929940, 0.5965649956210169);
          break;
        }
      }
    }; // class ShunnHamDriver<Simplex<3>>

  } // namespace Cubature
} // namespace FEAT

#endif // KERNEL_CUBATURE_SHUNN_HAM_DRIVER_HPP
